##合并临床资料



library(ConsensusClusterPlus)

d<- read.csv("Proportion_cluster8.csv", header = T,row.names = 1, check.names = F)



library(ClassDiscovery) #

Sys.setenv(LANGUAGE = "en") #显示英文报错信息
options(stringsAsFactors = FALSE) #禁止chr转成factor
mat <- d[,c(1,2,4)]
mat<- scale(  t(mat) )

mat[mat > 1] = 1
mat[mat < -1] = -1
mat[1:3,1:3]


annCol <-d[,c(11:18,27:30)]
annCol$GNAQ<-ifelse(annCol$GNAQ==1, "mutant","wildtype")
annCol$GNA11 <-ifelse(annCol$GNA11 ==1, "mutant","wildtype")
annCol$EIF1AX  <-ifelse(annCol$EIF1AX  ==1, "mutant","wildtype")
annCol$SF3B1  <-ifelse(annCol$SF3B1  ==1, "mutant","wildtype")
head(annCol)
colnames(annCol)
annCol[is.na(annCol) | annCol == ""] <- "N/A"
#annCol<-annCol[,3:10]
annColors <- list()
annColors[["gender"]] <- c("male"="#79B789", "female" = "#B5262A")
annColors[["vital_status"]] <- c("alive" = "#79B789", "dead" = "#B5262A")
annColors[["age"]] <- colorRampPalette(c("#00FFFF", "#2F4F4F"))(70)
annColors[["stage"]] <- c("Stage II"="#1E90FF","Stage III"="#C71585","Stage IV"="#DAA520")
annColors[["histological_type"]] <- c("epithelioid"="#1E90FF","spindle"="#C71585","mixed"="#DAA520")
annColors[["MetStatus"]] <- c("Metastatic"="#FFFF66","Non-metastatic"="#9C8C2E")
annColors[["chromosome.3.status"]] <- c("disomy"="blue","monosomy"="red")
annColors[["SCNA_Cluster"]]<-c("A"="#D848A9","B"="#D9681F","C"="#8BE36C","D"="#4169E1")
annColors[["GNAQ"]] <- c("mutant"="#000000","wildtype"="#DCDCDC")
annColors[["GNA11"]] <- c("mutant"="#000000","wildtype"="#DCDCDC")
annColors[["EIF1AX"]] <- c("mutant"="#000000","wildtype"="#DCDCDC")
annColors[["SF3B1"]] <- c("mutant"="#000000","wildtype"="#DCDCDC")


library(pheatmap)
library(RColorBrewer)
library(gplots)
pheatmap(mat,
         
         cluster_rows = F,
         cluster_cols = F,
         color = colorRampPalette(c("#343493", "white", "#C24A45"))(64),
         annotation_col = annCol,
         annotation_colors = annColors,
         show_rownames = T, show_colnames = F
)



hcs <- hclust(distanceMatrix(as.matrix(mat), "pearson"), "ward.D") # 请阅读distanceMatrix()以及hclust()，了解更多distance测度和linkage方法
hcg <- hclust(distanceMatrix(t(as.matrix(mat)), "pearson"), "ward.D") # 注意距离函数是针对列的，所以对行聚类要转置
group <- cutree(hcs,k=3)

# 增加一行annotation，及其配色
annCol$Clust2 <- paste0("Subtype",group[rownames(annCol)])
annColors[["Clust2"]] <- c("Subtype1"="#008ECB","Subtype2"="#EA921D","Subtype3"="#D14039")#,"Cluster4"="green")

pheatmap(mat,
         
         color =colorRampPalette(rev(brewer.pal(n = 8, name = "RdBu")))(100),
         cluster_rows = F,
         cluster_cols = hcs,
         annotation_col = annCol,
         annotation_colors = annColors,
         show_rownames = T,show_colnames = F
)

#把聚类后的样本顺序保存到文件
sample_order <- data.frame(row.names = seq(1:length(hcs$labels)), sample = hcs$labels, group = group)
sample_order <- sample_order[hcs$order,] #按聚类后的顺序排
sample_order$ori.order <- row.names(sample_order) #把最初的顺序保存在ori.order列里
write.csv(sample_order, "sample_order_TCGA.csv", quote = F, row.names = F)
head(annCol)



library(survival)
library(survminer)
#clinical<-read.csv("clinical_all.csv",row.names = 1,check.names = F,stringsAsFactors = F)

#clinical<-read.table("UVM_clinicalMatrix.txt",sep="\t",row.names=1,header = TRUE)
#clinical<-clinical[,c('OS.time','OS','vital_status')]
#write.csv(clinical, "clinic.csv", quote = F, row.names = T)


Sur <- Surv(d[,9],d[,10])
sfit <- survfit(Sur ~ annCol[,13],data=as.data.frame( annCol[,13]) )  
ggsurvplot(sfit, 
           conf.int=F, #置信区间
           #fun="pct",
           pval=TRUE,
           palette = c("#008ECB", "#EA921D","#D14039"),
           pval.method = T,
           risk.table =T, 
           ncensor.plot = F,
           surv.median.line="hv",
           legend.labs=c('Subtype1','Subtype2',"Subtype3"))+
  labs(x = "Days")


###HX_cohort
d<- read.csv("Proportion_HX.csv", header = T,row.names = 1, check.names = F)



mat <- d[,c(4,5,6)]
mat<- scale(  t(mat) )

mat[mat > 1] = 1
mat[mat < -1] = -1
mat[1:3,1:3]

d$vital_status<- ifelse(d$fustat==1,  'dead','alive')
annCol <-d[,c(11:13,15:20)]
annCol$MetStatus<- ifelse(annCol$MetStatus==1,  'Metastatic','Non-metastatic')
#annCol$vital_status<- ifelse(annCol$fustat==1,  'dead','alive')

annCol$GNAQ<-ifelse(annCol$GNAQ==1, "mutant","wildtype")
annCol$GNA11 <-ifelse(annCol$GNA11 ==1, "mutant","wildtype")
annCol$EIF1AX  <-ifelse(annCol$EIF1AX  ==1, "mutant","wildtype")
annCol$SF3B1  <-ifelse(annCol$SF3B1  ==1, "mutant","wildtype")
head(annCol)
colnames(annCol)
annCol[is.na(annCol) | annCol == ""] <- "N/A"
#annCol<-annCol[,3:10]
annColors <- list()
annColors[["gender"]] <- c("male"="#79B789", "female" = "#B5262A")
annColors[["vital_status"]] <- c("alive" = "#79B789", "dead" = "#B5262A")
annColors[["age"]] <- colorRampPalette(c("#00FFFF", "#2F4F4F"))(70)
#annColors[["stage"]] <- c("Stage II"="#1E90FF","Stage III"="#C71585","Stage IV"="#DAA520")
annColors[["histological_type"]] <- c("epithelioid"="#1E90FF","spindle"="#C71585","mix"="#DAA520")
annColors[["MetStatus"]] <- c("Metastatic"="#FFFF66","Non-metastatic"="#9C8C2E")
#annColors[["chromosome.3.status"]] <- c("disomy"="blue","monosomy"="red")
#annColors[["SCNA_Cluster"]]<-c("A"="#D848A9","B"="#D9681F","C"="#8BE36C","D"="#4169E1")
annColors[["GNAQ"]] <- c("mutant"="#000000","wildtype"="#DCDCDC")
annColors[["GNA11"]] <- c("mutant"="#000000","wildtype"="#DCDCDC")
annColors[["EIF1AX"]] <- c("mutant"="#000000","wildtype"="#DCDCDC")
annColors[["SF3B1"]] <- c("mutant"="#000000","wildtype"="#DCDCDC")


library(pheatmap)
library(gplots) 
pheatmap(mat,
         
         cluster_rows = F,
         cluster_cols = F,
         color = colorRampPalette(c("#343493", "white", "#C24A45"))(64),
         annotation_col = annCol,
         annotation_colors = annColors,
         show_rownames = T, show_colnames = F
)



hcs <- hclust(distanceMatrix(as.matrix(mat), "pearson"), "ward.D") # 请阅读distanceMatrix()以及hclust()，了解更多distance测度和linkage方法
hcg <- hclust(distanceMatrix(t(as.matrix(mat)), "pearson"), "ward.D") # 注意距离函数是针对列的，所以对行聚类要转置
group <- cutree(hcs,k=3)

# 增加一行annotation，及其配色
annCol$Clust2 <- paste0("Subtype",group[rownames(annCol)])
annColors[["Clust2"]] <- c("Subtype1"="#D14039","Subtype2"="#008ECB","Subtype3"="#EA921D")#,"Cluster4"="green")

pheatmap(mat,
         
         color =colorRampPalette(rev(brewer.pal(n = 8, name = "RdBu")))(100),
         cluster_rows = F,
         cluster_cols = hcs,
         annotation_col = annCol,
         annotation_colors = annColors,
         show_rownames = T,show_colnames = F
)

#把聚类后的样本顺序保存到文件
sample_order <- data.frame(row.names = seq(1:length(hcs$labels)), sample = hcs$labels, group = group)
sample_order <- sample_order[hcs$order,] #按聚类后的顺序排
sample_order$ori.order <- row.names(sample_order) #把最初的顺序保存在ori.order列里
write.csv(sample_order, "sample_order_HX.csv", quote = F, row.names = F)
head(annCol)



library(survival)
library(survminer)
#clinical<-read.csv("clinical_all.csv",row.names = 1,check.names = F,stringsAsFactors = F)

#clinical<-read.table("UVM_clinicalMatrix.txt",sep="\t",row.names=1,header = TRUE)
#clinical<-clinical[,c('OS.time','OS','vital_status')]
#write.csv(clinical, "clinic.csv", quote = F, row.names = T)


Sur <- Surv(d[,9]*30,d[,10])
sfit <- survfit(Sur ~ annCol[,10],data=as.data.frame( annCol[,10]) )  
ggsurvplot(sfit, 
           conf.int=F, #置信区间
           #fun="pct",
           pval=TRUE,
           palette = c("#008ECB", "#D14039","#EA921D"),
           pval.method = T,
           risk.table =T, 
           ncensor.plot = F,
           surv.median.line="hv",
           legend.labs=c('Subtype1','Subtype2',"Subtype3"))+
  labs(x = "Days")




###TCGA_cohort
LIRI<- read.csv("Proportion_cluster8.csv", header = T,row.names = 1, check.names = F)
surv.dat<-LIRI[,c(1,2,4,9,10)]

f = survival::coxph(Surv(futime, fustat) ~ ., data = surv.dat)
#  2.risk point nomgram.points and lp
riskscore = f$linear.predictors




##TCGA

library(survival)
library(survivalROC)
library(survminer)



riskscore <- data.frame(riskscore = as.numeric(riskscore),
                        #group = ifelse(riskscore >1.2,"HRisk","LRisk"), # 根据中位数分组
                        row.names = rownames(surv.dat),
                        OS.time = surv.dat$futime,
                        OS = surv.dat$fustat,
                        stringsAsFactors = F)

res.cut <- surv_cutpoint(riskscore, time = "OS.time",
                         event = "OS",
                         variables = names(riskscore)[1],
                         minprop = 0.3)

res.cat <- surv_categorize(res.cut)
riskscore$group<-res.cat$riskscore

Sur <- Surv(riskscore[,2] ,riskscore[,3])
sfit <- survfit(Sur ~ riskscore[,4],data=as.data.frame( riskscore[,4]) )
ggsurvplot(sfit,
           conf.int=F, #置信区间
           #fun="pct",
           pval=TRUE,
           palette = c("#EA921D", "#008ECB"),
           pval.method = T,
           risk.table =T,
           ncensor.plot = F,
           surv.median.line="hv",
           legend.labs=c('High_risk','Low_risk'))+
  labs(x = "Days")


###HX_cohort
LIRI<- read.csv("Proportion_HX.csv", header = T,row.names = 1, check.names = F)

surv.dat<-LIRI[,c(4:6,9,10)]

f = survival::coxph(Surv(futime, fustat) ~ ., data = surv.dat)
#  2.risk point nomgram.points and lp
riskscore = f$linear.predictors



library(survival)
library(survivalROC)
library(survminer)



riskscore <- data.frame(riskscore = as.numeric(riskscore),
                        #group = ifelse(riskscore >1.2,"HRisk","LRisk"), # 根据中位数分组
                        row.names = rownames(surv.dat),
                        OS.time = surv.dat$futime,
                        OS = surv.dat$fustat,
                        stringsAsFactors = F)

res.cut <- surv_cutpoint(riskscore, time = "OS.time",
                         event = "OS",
                         variables = names(riskscore)[1],
                         minprop = 0.3)

res.cat <- surv_categorize(res.cut)
riskscore$group<-res.cat$riskscore

Sur <- Surv(riskscore[,2]*30 ,riskscore[,3])
sfit <- survfit(Sur ~ riskscore[,4],data=as.data.frame( riskscore[,4]) )
ggsurvplot(sfit,
           conf.int=F, #置信区间
           #fun="pct",
           pval=TRUE,
           palette = c("#EA921D", "#008ECB"),
           pval.method = T,
           risk.table =T,
           ncensor.plot = F,
           surv.median.line="hv",
           legend.labs=c('High_risk','Low_risk'))+
  labs(x = "Days")
